//-------------------------------------------------------------------------//
//                                                                         //
//  This benchmark is an OpenMP C version of the NPB LU code. This OpenMP  //
//  C version is developed by the Center for Manycore Programming at Seoul //
//  National University and derived from the OpenMP Fortran versions in    //
//  "NPB3.3-OMP" developed by NAS.                                         //
//                                                                         //
//  Permission to use, copy, distribute and modify this software for any   //
//  purpose with or without fee is hereby granted. This software is        //
//  provided "as is" without express or implied warranty.                  //
//                                                                         //
//  Information on NPB 3.3, including the technical report, the original   //
//  specifications, source code, results and information on how to submit  //
//  new results, is available at:                                          //
//                                                                         //
//           http://www.nas.nasa.gov/Software/NPB/                         //
//                                                                         //
//  Send comments or suggestions for this OpenMP C version to              //
//  cmp@aces.snu.ac.kr                                                     //
//                                                                         //
//          Center for Manycore Programming                                //
//          School of Computer Science and Engineering                     //
//          Seoul National University                                      //
//          Seoul 151-744, Korea                                           //
//                                                                         //
//          E-mail:  cmp@aces.snu.ac.kr                                    //
//                                                                         //
//-------------------------------------------------------------------------//

//-------------------------------------------------------------------------//
// Authors: Sangmin Seo, Jungwon Kim, Jun Lee, Jeongho Nah, Gangwon Jo,    //
//          and Jaejin Lee                                                 //
//-------------------------------------------------------------------------//

#include <stdio.h>
#include "applu.incl"
#include <stdlib.h>

#include <mpi.h>

void pintgr()
{



  //---------------------------------------------------------------------
  // local variables
  //---------------------------------------------------------------------
  int i, j, k;
  int ibeg, ifin, ifin1;
  int jbeg, jfin, jfin1;
  double phi1[ISIZ3+2][ISIZ2+2];
  double phi2[ISIZ3+2][ISIZ2+2];
  double frc1, frc2, frc3;

  //---------------------------------------------------------------------
  // set up the sub-domains for integeration in each processor
  //---------------------------------------------------------------------
  ibeg = ii1;
  ifin = ii2;
  jbeg = ji1;
  jfin = ji2;
  ifin1 = ifin - 1;
  jfin1 = jfin - 1;


MPI_Status status;


   int my_id = 0, id_num = 0,upper_bound = 0, lower_bound = 0, pros_num = 0, average_span = 0, root_process=0;
   int ierr, upper_bound_two = 0, lower_bound_two = 0, average_span_two = 0,upper_bound_three = 0, lower_bound_three = 0, average_span_three = 0;

   MPI_Comm comm;

   ierr = MPI_Comm_rank(MPI_COMM_WORLD, &my_id);
   ierr = MPI_Comm_size(MPI_COMM_WORLD, &pros_num);
   
   int *delta = (int *) malloc(pros_num*sizeof(int));
   int *disp = (int *) malloc(pros_num*sizeof(int));

   double scount = 0;
   int rcount[pros_num];
   int countit;
   


   root_process = 0;

   MPI_Request request, request2;
  

   average_span = (jfin-jbeg)/pros_num;
   average_span_two = (jfin1-jbeg)/pros_num;
   average_span_three = ((ki2-1)-ki1)/pros_num;   

        if (my_id == root_process){

        lower_bound = jbeg;
        upper_bound = average_span;

        lower_bound_two = jbeg;
        upper_bound_two = average_span_two;

        lower_bound_three = ki1;
        upper_bound_three = average_span_three+1;


           if ((my_id+1)  == pros_num){
                if(upper_bound < jfin-1 ){
                        upper_bound = jfin-1;
                   }
                if(upper_bound_two < jfin1-1 ){
                        upper_bound_two = jfin1-1;
                   }
                if(upper_bound_three < ki2-2 ){
                        upper_bound_three = ki2-2;
                   }

                }

         }
        else{

        lower_bound = (average_span *  my_id) + 1;
        upper_bound = average_span * (my_id+1);

        lower_bound_two = (average_span_two *  my_id) + 1;
        upper_bound_two = average_span_two * (my_id+1);

        lower_bound_three = (average_span_three *  my_id) + 2;
        upper_bound_three = (average_span_three * (my_id+1))+1;

          if ((my_id+1) == pros_num){
                if(upper_bound < jfin-1 ){
                        upper_bound = jfin-1;
                   }
                if(upper_bound_two < jfin1-1 ){
                        upper_bound_two = jfin1-1;
                   }
                if(upper_bound_three < ki2-2 ){
                        upper_bound_three = ki2-2;
                   }

                 }
             }





//  #pragma omp parallel default(shared) private(i,j,k) \
                       shared(ki1,ki2,ifin,ibeg,jfin,jbeg,ifin1,jfin1)
  {

//printf(" j is from  %u to %u and i is from %u to %u and lower is is %u and upper is %u\n",jbeg,jfin,ibeg,ifin,lower_bound,upper_bound);
//  #pragma omp for nowait
  for (j = jbeg; j < jfin; j++) {

//for (j = lower_bound; j <= upper_bound; j++) {
    for (i = ibeg; i < ifin; i++) {
      k = ki1;

      phi1[j][i] = C2*(  u[k][j][i][4]
          - 0.50 * (  u[k][j][i][1] * u[k][j][i][1]
                    + u[k][j][i][2] * u[k][j][i][2]
                    + u[k][j][i][3] * u[k][j][i][3] )
                   / u[k][j][i][0] );

      k = ki2 - 1;

      phi2[j][i] = C2*(  u[k][j][i][4]
          - 0.50 * (  u[k][j][i][1] * u[k][j][i][1]
                    + u[k][j][i][2] * u[k][j][i][2]
                    + u[k][j][i][3] * u[k][j][i][3] )
                   / u[k][j][i][0] );
    }
  }




/*

for (countit=0;countit < pros_num; countit++){

//disp[countit]=0;
rcount[countit] = delta[countit]*(ISIZ2+2);
disp[countit] = disp[countit]*(ISIZ2+2);
//rcount[countit]=scount;
printf("\n\nrcount is %d and disp is %d\n\n",rcount[countit],disp[countit]);
 }

//double scount = (ISIZ3+2)*(ISIZ2+2);
//double rcount = scount_1;
//printf("......................i get here at least");

ierr = MPI_Allgatherv(MPI_IN_PLACE,1,MPI_DOUBLE,&phi1[1][1],rcount,disp,MPI_DOUBLE,MPI_COMM_WORLD);
ierr = MPI_Allgatherv(MPI_IN_PLACE,1,MPI_DOUBLE,&phi2[1][1],rcount,disp,MPI_DOUBLE,MPI_COMM_WORLD);

//ierr=MPI_Allgatherv( &phi1[lower_bound][1], scount, MPI_DOUBLE, &phi1[1][1], rcount,disp, MPI_DOUBLE, MPI_COMM_WORLD);
//ierr=MPI_Allgatherv( &phi2[lower_bound][1], scount, MPI_DOUBLE, &phi2[1][1], rcount,disp, MPI_DOUBLE, MPI_COMM_WORLD);

//double scount_two =scount;// (average_span_two)*(ISIZ2+2);
//double rcount_two = scount;


if(my_id == 0){
printf("\n ..........**********<<<<<<<<for root upper bound is %u and  lower bound %u and scount is %f\n ",upper_bound, lower_bound, scount );
//ierr=MPI_Allgather( &phi1[lower_bound][1], scount, MPI_DOUBLE, &phi1[1][1], scount, MPI_DOUBLE, MPI_COMM_WORLD);
//ierr=MPI_Allgather( &phi2[lower_bound][1], scount, MPI_DOUBLE, &phi2[1][1], scount, MPI_DOUBLE, MPI_COMM_WORLD);
ierr = MPI_Isend( &phi1[lower_bound][1], scount , MPI_DOUBLE,1, MPI_ANY_TAG, MPI_COMM_WORLD, &request);
ierr = MPI_Irecv(&phi1[lower_bound][1],scount_1 , MPI_DOUBLE,1, MPI_ANY_TAG, MPI_COMM_WORLD, &request);    

ierr = MPI_Isend( &phi2[lower_bound][1], scount , MPI_DOUBLE,1, MPI_ANY_TAG, MPI_COMM_WORLD, &request);
ierr = MPI_Irecv(&phi2[lower_bound][1],scount_1 , MPI_DOUBLE,1, MPI_ANY_TAG, MPI_COMM_WORLD, &request);

}

if(my_id == 1){
printf("\n..........**********<<<<<<<< for non root upper bound is %u and  lower bound %u and scount is %f\n ",upper_bound, lower_bound, scount );
//ierr=MPI_Allgather( &phi1[lower_bound][1], scount_1, MPI_DOUBLE, &phi1[1][1], scount_1, MPI_DOUBLE, MPI_COMM_WORLD);
//ierr=MPI_Allgather( &phi2[lower_bound][1], scount_1, MPI_DOUBLE, &phi2[1][1], scount_1, MPI_DOUBLE, MPI_COMM_WORLD);
//ierr = MPI_Recv( &phi1[lower_bound][1], num_rows_to_receive, MPI_INT,root_process, send_data_tag, MPI_COMM_WORLD, &status);
ierr = MPI_Isend( &phi1[lower_bound][1], scount_1 , MPI_DOUBLE,0, MPI_ANY_TAG, MPI_COMM_WORLD, &request);
ierr = MPI_Irecv(&phi1[lower_bound][1],scount , MPI_DOUBLE,0, MPI_ANY_TAG, MPI_COMM_WORLD, &request);

ierr = MPI_Isend( &phi2[lower_bound][1], scount_1 , MPI_DOUBLE,0, MPI_ANY_TAG, MPI_COMM_WORLD, &request);
ierr = MPI_Irecv(&phi2[lower_bound][1],scount , MPI_DOUBLE,0, MPI_ANY_TAG, MPI_COMM_WORLD, &request);
}


//MPI_Bcast( phi1, scount, MPI_DOUBLE, root_process, MPI_COMM_WORLD); 
//MPI_Bcast( phi2, scount_two, MPI_DOUBLE, root_process, MPI_COMM_WORLD);
//....con here..^

int MPI_Barrier( MPI_Comm comm );
*/

//  #pragma omp single
if(my_id == root_process){
  frc1 = 0.0;
}

MPI_Bcast( &frc1, 0.0, MPI_DOUBLE, root_process, MPI_COMM_WORLD); 



//  #pragma omp for reduction(+:frc1)
//  for (j = jbeg; j < jfin1; j++) {

  for (j = lower_bound_two; j <= upper_bound_two; j++) {
    for (i = ibeg; i < ifin1; i++) {
      frc1 = frc1 + (  phi1[j][i]
                     + phi1[j][i+1]
                     + phi1[j+1][i]
                     + phi1[j+1][i+1]
                     + phi2[j][i]
                     + phi2[j][i+1]
                     + phi2[j+1][i]
                     + phi2[j+1][i+1] );
    }
  }


//ierr=MPI_Allreduce( &frc1,&frc1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
ierr=MPI_Allreduce( MPI_IN_PLACE,&frc1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

//  #pragma omp single nowait
if(my_id == root_process){
  frc1 = dxi * deta * frc1;
}
MPI_Bcast( &frc1, 1, MPI_DOUBLE, root_process, MPI_COMM_WORLD);
//int MPI_Barrier( MPI_Comm comm );


//  #pragma omp for nowait
  for (k = ki1; k < ki2; k++) {
    for (i = ibeg; i < ifin; i++) {
      phi1[k][i] = C2*(  u[k][jbeg][i][4]
          - 0.50 * (  u[k][jbeg][i][1] * u[k][jbeg][i][1]
                    + u[k][jbeg][i][2] * u[k][jbeg][i][2]
                    + u[k][jbeg][i][3] * u[k][jbeg][i][3] )
                   / u[k][jbeg][i][0] );
    }
  }

//  #pragma omp for nowait
  for (k = ki1; k < ki2; k++) {
    for (i = ibeg; i < ifin; i++) {
      phi2[k][i] = C2*(  u[k][jfin-1][i][4]
          - 0.50 * (  u[k][jfin-1][i][1] * u[k][jfin-1][i][1]
                    + u[k][jfin-1][i][2] * u[k][jfin-1][i][2]
                    + u[k][jfin-1][i][3] * u[k][jfin-1][i][3] )
                   / u[k][jfin-1][i][0] );
    }
  }

//  #pragma omp single
  frc2 = 0.0;

//printf("\nki1 is %u to ki2-1 is %u and cal is from %u to %u.....>>%u\n",ki1,ki2-1,lower_bound_three,upper_bound_three,my_id);
//  #pragma omp for reduction(+:frc2)
//  for (k = ki1; k < ki2-1; k++) {

for (k = lower_bound_three; k <= upper_bound_three; k++) {

    for (i = ibeg; i < ifin1; i++) {
      frc2 = frc2 + (  phi1[k][i]
                     + phi1[k][i+1]
                     + phi1[k+1][i]
                     + phi1[k+1][i+1]
                     + phi2[k][i]
                     + phi2[k][i+1]
                     + phi2[k+1][i]
                     + phi2[k+1][i+1] );
    }
  }

//ierr=MPI_Allreduce( &frc2,&frc2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
ierr=MPI_Allreduce( MPI_IN_PLACE,&frc2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

if(my_id == root_process){
  frc2 = dxi * dzeta * frc2;
}
MPI_Bcast( &frc2, 1, MPI_DOUBLE, root_process, MPI_COMM_WORLD);

//  #pragma omp single nowait
//  frc2 = dxi * dzeta * frc2;

//  #pragma omp for nowait
  for (k = ki1; k < ki2; k++) {
    for (j = jbeg; j < jfin; j++) {
      phi1[k][j] = C2*(  u[k][j][ibeg][4]
          - 0.50 * (  u[k][j][ibeg][1] * u[k][j][ibeg][1]
                    + u[k][j][ibeg][2] * u[k][j][ibeg][2]
                    + u[k][j][ibeg][3] * u[k][j][ibeg][3] )
                   / u[k][j][ibeg][0] );
    }
  }

//  #pragma omp for nowait
  for (k = ki1; k < ki2; k++) {
    for (j = jbeg; j < jfin; j++) {
      phi2[k][j] = C2*(  u[k][j][ifin-1][4]
          - 0.50 * (  u[k][j][ifin-1][1] * u[k][j][ifin-1][1]
                    + u[k][j][ifin-1][2] * u[k][j][ifin-1][2]
                    + u[k][j][ifin-1][3] * u[k][j][ifin-1][3] )
                   / u[k][j][ifin-1][0] );
    }
  }

//  #pragma omp single
  frc3 = 0.0;

//  #pragma omp for reduction(+:frc3)
//  for (k = ki1; k < ki2-1; k++) {

for (k = lower_bound_three; k <= upper_bound_three; k++) {
    for (j = jbeg; j < jfin1; j++) {
      frc3 = frc3 + (  phi1[k][j]
                     + phi1[k][j+1]
                     + phi1[k+1][j]
                     + phi1[k+1][j+1]
                     + phi2[k][j]
                     + phi2[k][j+1]
                     + phi2[k+1][j]
                     + phi2[k+1][j+1] );
    }
  }
//ierr=MPI_Allreduce( &frc3,&frc3, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
ierr=MPI_Allreduce(MPI_IN_PLACE,&frc3, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

if(my_id == root_process){
frc3 = deta * dzeta * frc3;
}

MPI_Bcast( &frc3, 1, MPI_DOUBLE, root_process, MPI_COMM_WORLD);
//  #pragma omp single nowait
//  frc3 = deta * dzeta * frc3;
  } //end parallel

  frc = 0.25 * ( frc1 + frc2 + frc3 );
  //printf("\n\n     surface integral = %12.5E\n\n\n", frc);
}

